#include "afdcpp.h"
#include "utils.h"
#include "renorm_dia.h"
#include "misc.h"
#include "feynrules.h"
// print particle content
void p_table()
{
 cout<<setw(17)<<"Field name"<<"field number"<<endl;
  cout<<setw(17)<<"No field"<<no_field<<endl;
  cout<<setw(17)<<"phi-scalar"<<phi<<endl;
  cout<<setw(17)<<"e-electron"<<e<<endl;
  cout<<setw(17)<<"E-positron"<<E<<endl;
  cout<<setw(17)<<"A-photon"<<A<<endl;
}

/// Non member template functions


///true if vertex_degree of graph less then
// vertex_degree of root graph
template <typename T>
bool is_local_ext_vert(T vdis,Graph& Gtmp)
{
  if(out_degree(vdis,Gtmp))
    {
      if(Gtmp.is_root())return false;
      else  return((out_degree(vdis,Gtmp) < out_degree( Gtmp.local_to_global(vdis), Gtmp.root())));
    }
  else return false;
}

/// replacement for Graph::find_vertex()
// synpsys:
//  vdis to find,global, local graph
template <typename T>
bool find_vert(T vdis,Graph& Gtmp)
{
  Graph::vertex_iterator vi,vi_end;
  for(tie(vi,vi_end) = vertices(Gtmp);vi != vi_end;++vi)
    if(Gtmp.local_to_global(*vi) == vdis)return true;
  return false;
}
/// analog of Graph::find_vertex()
// synpsys:
//  edis to find,global, local graph
/*template <typename T>
bool find_edge(T edis,Graph& Gtmp)
{
  Graph::edge_iterator ei,ei_end;
  for(tie(ei,ei_end) = edges(Gtmp);ei != ei_end;++ei)
    if(Gtmp.local_to_global(*ei) == edis)return true;
  return false;
}*/
/// Check if vertex is at one of ends of spinor line
template <typename T>
bool is_ext_spin(T vdis,Graph& Gtmp)
{
  Graph::out_edge_iterator oei,oei_end;
  int n_sp = 0;
  for(tie(oei,oei_end) = out_edges(vertex(vdis,Gtmp),Gtmp);oei!=oei_end;++oei)
    if(Gtmp[*oei].type == 1) n_sp++;
  cout<<" n_sp++ "<<n_sp<<endl;
  return (is_external_vert(vdis,Gtmp) && (1 == n_sp) );    
}


/*
/// Get information about field of GLOBAL external vertex
template <typename T>
std::pair<field,bool> get_global_ext_field_info(Graph::vertex_descriptor,Graph& G)
{
}
*/





/// From list like G1,G2,G3...Gn
/// creates list of terms like G1*G2,G1*G3...
 /*ex afdcpp::get_terms(Graphref_lst l1)  
{
  std::vector<bool> tmk(l1.size(),false);
  typedef std::list<Graph> M_lst;
  std::list<M_lst> lMl;
  ex o_ex;
  //  inc(tmk);
  while(!is_full(tmk)) 
    {
      inc(tmk);
      M_lst ml;
      int cnt = 0;
      BOOST_FOREACH(Graph g,l1)
        {
          if(tmk[cnt])ml.push_back(g);
          cnt++;
        }
      cout<<"Extracting terms "<<endl;
      ex cvb = term_ex(ml);
      o_ex+=cvb;
      cout<<colors::cyan<<"final expr: "<<cvb<<colors::def_col<<endl;
      lMl.push_back(ml);
      //   inc(tmk);
    }

  cout<<" number of permutations(lMl size) "<<lMl.size()<<endl;
  cout<<colors::red<<"Expr for substraction:"<<colors::def_col<<endl;
  cout<<o_ex<<endl;
  return o_ex;
}


 */


 
// <div_graph, T(div_graph)>
std::pair<Graph,Graph> afdcpp::K2(Graph& Gdiv,Graph& GT)
{
  /*  cout<<colors::blue<<"test in K"<<endl;
  print(GT);
  cout<<colors::def_col;

  Graph_1PI_lst l_div = make_subgraphs(GT);
  /*
  BOOST_FOREACH(Graph g,l_div)
    {
      if(isSubgraph(Gdiv,g))
        {
          cout<<"45 found div sub"<<endl;
          return g;
        }
    }
  *//*
  if(l_div.size())
    {
      for(Graph_1PI_lst::iterator it = l_div.begin();it!=l_div.end();++it)
        if(isSubgraph(Gdiv,*it))
          {
            cout<<"K2 found div sub"<<endl;
            //Graph g76 = *it;
            print_graph(*it);
            Graph tmtmgg = *it;///erefe
            return make_pair(tmtmgg,T(*it));
          }
    }
  else
    {
      cout<<" MEGAERROR=)"<<endl;
      return make_pair(Gdiv,GT);
    }
*/
}

exvector afdcpp::loop_momentum(ex dia_ex,exvector em_vec)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | Find loop momentums in dia expr     | ///
/// | and in external momentums vector    | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  exvector loop_lst;
  loop_lst.push_back(k1);
  loop_lst.push_back(k2);
  loop_lst.push_back(k3);
  loop_lst.push_back(k4);
  loop_lst.push_back(k5);
  loop_lst.push_back(k6);
  exvector ext_mom;
  std::list<ex> dia_mom;
  BOOST_FOREACH(ex mom,loop_lst)
    {
      if(dia_ex.has(mom))
        {
          cout<<"dia has "<<mom<<endl;
          dia_mom.push_back(mom);
        }
    }
  BOOST_FOREACH(ex external_mom,em_vec)
    {
      BOOST_FOREACH(ex mom,loop_lst)
        {
          if(external_mom.has(mom))
            {
              cout<<"ext has "<<mom<<endl;
              dia_mom.remove(mom);
            }          
        }
    }
  BOOST_FOREACH(ex m23,dia_mom)
    {
      cout<<"LOOP MOMENTUM: "<<m23<<endl;
    }
  exvector out_vec(dia_mom.begin(),dia_mom.end());
  return out_vec;
}


ex afdcpp::int_graph(Graph G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |Integral of graph in terms of epsilon| ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  ex out_e = ex_from_graph(G);
  exvector ev1 = get_ext_momentum(G);
  // One-loop integral function
  symbol rho("rho");
  ex C = Scalar2Pt(p1, mS, mS, 1, 1, rho);//, eps);
  // Perform epsilon-expansion
  C = C.series(eps == 0, 1);
  // Extract the divergent and finite parts
  ex fex = out_e.numer()*C.subs(lst(rho == 0,mS==0,p1==1)).normal();

  return fex;
}

ex afdcpp::div_part(Graph& G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | Divergent part of integral          | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  // ex projected_expr = projector(G);
  return  0;//projector(G);
}


ex afdcpp::term_ex(std::list<Graph>& l_in)
{
  cout<<"extraction procedure started for l_in = "<<l_in.size()<<" terms"<<endl;
  if(l_in.size() > 1)
    {
      ///searchig for zero terms
      for (std::list<Graph>::iterator it = l_in.begin();it != l_in.end();++it)
        for (std::list<Graph>::iterator it2 = l_in.begin();it2 != l_in.end();++it2)
          if((it != it2)&&(isStrongIntersected(*it,*it2))) 
            {
              cout<<"ZERO TERM"<<endl;
              return 0;
            }
      /// sorting the list
      cout<<"Not sorted";
      BOOST_FOREACH(Graph g,l_in)
        {
          cout<<" "<<num_vertices(g)<<" ";
        }
      cout<<endl;
      l_in.sort(isSubgraph);
      //      l_in.reverse();
      cout<<"After sorting";
      BOOST_FOREACH(Graph g3,l_in)
        {
          cout<<" "<<num_vertices(g3)<<" ";
        }
      cout<<endl;
      std::list<Graph> new_list;
      Graph G_tmp;
      Graph g_t;
      std::list<Graph> gl_to_store;
      for(std::list<Graph>::reverse_iterator rit = l_in.rbegin();rit != l_in.rend();++rit)
        {
          if(rit == l_in.rbegin()) G_tmp = (*rit);
          else if(rit == (--l_in.rend())) 
            {
              cout<<colors::red<<"T of last graph"<<endl;
              cout<<colors::def_col;
              gl_to_store.push_back(T(G_tmp));
              tie(G_tmp,g_t) = K2((*rit),gl_to_store.back());
              cout<<"gtmptest"<<endl;
              print_graph(G_tmp);
            }
          else
            {
              cout<<"Subgraphing "<<get_property(T(new_list.back()),graph_num)<<endl;
              cout<<endl<<endl<<"     T for K"<<endl<<endl<<endl;
              //              print(T(new_list.back()));
              gl_to_store.push_back(T(G_tmp));
              G_tmp = K((*rit),gl_to_store.back());
            }
          new_list.push_back(G_tmp);
        }
      cout<<"nmt : "<<new_list.size()<<endl;
      //      BOOST_FOREACH(Graph g7,new_list)
      for(std::list<Graph>::iterator it = new_list.begin();it!=new_list.end();++it)
        {
          cout<<"Graph:"<<endl;
          cout<<endl<<colors::red<<"   "<<endl;//<<
          cout<<ex_from_graph(*it)<<endl<<endl;
          print_graph(*it);
          if(!(*it).is_root())
            {
              cout<<"It's root"<<endl;
              print_graph((*it).root());
            }
        }
      cout<<colors::def_col;
      ex out_ex = ex_from_graph(g_t);
      return out_ex;
      //return 0;
    }
  else
    {
      cout<<"Working with single term"<<endl;
      // return ex_from_graph(*(l_in.begin()));
      return div_part(*(l_in.begin()));//*int_graph(T(*(l_in.begin())));
    }
}

/// Constructing program Functions 
scalar_products afdcpp::global_sp()
{
  scalar_products sp;
  sp.add(k1, p1, k1_p1);
  sp.add(k1, p2, k1_p2);
  sp.add(k1, q1, k1_q1);
  sp.add(k1, q2, k1_q2);
  sp.add(k2, p1, k2_p1);
  sp.add(k2, p2, k2_p2);
  sp.add(k2, q1, k2_q1);
  sp.add(k2, q2, k2_q2);
  sp.add(k3, p1, k3_p1);
  sp.add(k3, p2, k3_p2);
  sp.add(k3, q1, k3_q1);
  sp.add(k3, q2, k3_q2);
  sp.add(k4, p1, k4_p1);
  sp.add(k4, p2, k4_p2);
  sp.add(k4, q1, k4_q1);
  sp.add(k4, q2, k4_q2);
  sp.add(k5, p1, k5_p1);
  sp.add(k5, p2, k5_p2);
  sp.add(k5, q1, k5_q1);
  sp.add(k5, q2, k5_q2);
  sp.add(k6, p1, k6_p1);
  sp.add(k6, p2, k6_p2);
  sp.add(k6, q1, k6_q1);
  sp.add(k6, q2, k6_q2);



  sp.add(k1, k1, k1_k1);
  sp.add(k1, k2, k1_k2);
  sp.add(k1, k3, k1_k3);
  sp.add(k1, k4, k1_k4);
  sp.add(k1, k5, k1_k5);
  sp.add(k1, k6, k1_k6);
  sp.add(k2, k2, k2_k2);
  sp.add(k2, k3, k2_k3);
  sp.add(k2, k4, k2_k4);
  sp.add(k2, k5, k2_k5);
  sp.add(k2, k6, k2_k6);
  sp.add(k3, k3, k3_k3);
  sp.add(k3, k4, k3_k4);
  sp.add(k3, k5, k3_k5);
  sp.add(k3, k6, k3_k6);
  sp.add(k4, k4, k4_k4);
  sp.add(k4, k5, k4_k5);
  sp.add(k4, k6, k4_k6);
  sp.add(k5, k5, k5_k5);
  sp.add(k5, k6, k5_k6);
  sp.add(k6, k6, k6_k6);

  sp.add(p1,p1, p1_p1);
  sp.add(p2,p2, p2_p2);
  sp.add(q1,q1, q1_q1);
  sp.add(q2,q2, q2_q2);
  sp.add(p1,q1, p1_q1);
  sp.add(p2,q2, p1_q2);
  sp.add(p2,q1, p2_q1);
  sp.add(p2,q2, p2_q2);

 

  //  sp.add(p1, p1, pow(m,2)); // A and B are orthogonal
  //  sp.add(k1, k1, 0); // A and C are orthogonal
  //  sp.add(k2, k2, 0); // A^2 = 4 (A has length 2)
  //  sp.add(k1, p1, k1_p1);
  // sp.add(k1,k2, k1_k2);
  //  sp.add(k2, p1, k2_p1);
  return sp;
}





/// DRAW DIAGRAM TO FILE

/*void afdcpp::draw_lst(Graph_1PI_lst& l)
{
  std::ofstream of("draw/out.tex");
  of<<"\\documentclass{article}\n\\usepackage{feynmf}%%%{feynmp}\n\\unitlength=1mm\n\\begin{document}\n";
  
  BOOST_FOREACH(Graph g,l)
    {
      of<<endl<<"\\\\"<<"\n    dia\n"<<"\\\\"<<endl;
      of<<draw(g)<<endl;
    }
  of<<"\\end{document}\n";
}
*/
void afdcpp::draw2(Graph& G,FD::page& pg)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | FeynDiagram drawing procedure       | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{


 Graph::edge_iterator ei,ei_end;
 for (tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) {

     put(edge_weight, G, *ei, 0.00001);
 }
 // bool ok = kamada_kawai_spring_layout(G,get(vertex_position, G),get(edge_weight, G),square_topology<>(1.0),edge_length(50.0),boost::layout_tolerance<>(10), 1);
                                      //  BOOST_CHECK(ok);

 //random_graph_layout
 //circle_graph_layout
    fruchterman_reingold_force_directed_layout
      (G,get(vertex_position, G),square_topology<>(20.0));
 using FD::page; 
 // page pg;
 FD::FeynDiagram fd(pg);

 std::map<Graph::vertex_descriptor,FD::vertex_dot*> vmap;
 Graph::vertex_iterator vi,vi_end;
 for (tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi) {

   cout<<"YYYYYYYYYYYYYY x,y=("<<get(vertex_position, G)[*vi][0]<<", "<<get(vertex_position, G)[*vi][1]<<endl;
   FD::vertex_dot* v1 = new FD::vertex_dot(fd,get(vertex_position, G)[*vi][0],get(vertex_position, G)[*vi][1]);
   vmap[*vi] = v1;
 }

 for (tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
   {
     if(G[*ei].type == 1) FD::line_plain* f1 = new FD::line_plain(fd,*vmap[source(*ei,G)],*vmap[target(*ei,G)]);
     else if(G[*ei].type == 0)
       {
         FD::line_wiggle* f1 = new FD::line_wiggle(fd,*vmap[source(*ei,G)],*vmap[target(*ei,G)]);
         FD::xy midpt, arcpt;
         midpt = 0.5 * (*vmap[source(*ei,G)] + *vmap[target(*ei,G)]);
         arcpt = midpt + 4.0 * (*vmap[source(*ei,G)] - *vmap[target(*ei,G)]).perp();
         Graph::edge_descriptor edis;
         //bool hp = has_par_edges(edis((*ei)),G);
         bool hp = true;
         if(hp) f1->arcthru(arcpt);
       }
     else if(G[*ei].type == -1)
       { FD::line_plain* f1 = new FD::line_plain(fd,*vmap[source(*ei,G)],*vmap[target(*ei,G)]);
         f1->dashon.settrue();
       }
     else cout<<"err detected"<<endl;
   }

 std::cout << "Triangle layout (Kamada-Kawai).\n";
 //print_graph_layout(g, get(vertex_position, g));
 
 


 // return pg;
}

std::string afdcpp::draw(Graph& G)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | Feynmf package drowing procedure    | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  if(!G.is_root()){
  //random_number_generator<boost::mt19937> rng; 
 
  std::ostringstream of;
  //write_graphviz(of, G);
  //  of<<"\\begin{fmffile}{fmftempl"<<get_property(G,graph_num)<<"}\n\\begin{fmfchar*}(80,50)\n";
  of<<"\\begin{fmffile}{fmftempl"<<draw_rand()<<"}\n\\begin{fmfchar*}(80,50)\n";

  // add out spinor lines and vertices with external legs 
  Graph::vertex_iterator vi,vi_end;
  for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi) G[*vi].visited = false;
  Graph::edge_iterator ei,ei_end;
  for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) G[*ei].visited = false;
  Graph::vertex_descriptor vdis;
  Graph::edge_descriptor edis;
  
  /// count external vertices
  std::list<Graph::vertex_descriptor> vert_lst;
  for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
    {
      Graph::vertex_descriptor vd = *vi;
    if(is_external_vert(vd,G)) vert_lst.push_back(*vi);
    }
  //  of<<"num ext: "<<ext_vert_lst.size();
  /// mark vertices as left and right
  int cntr = 0;
  std::string str_left,str_right,left_str_tag;
  bool h1 = false,h2 = false;
  for(std::list<Graph::vertex_descriptor>::iterator it = vert_lst.begin(); it != vert_lst.end();++it)
    {
          
      cntr++;
      if(cntr <= vert_lst.size()/2) 
        if(cntr == 1) str_left+=lexical_cast<string>(*it);
        else str_left+=","+lexical_cast<string>(*it);
      else
        if(!h2){
          str_right+=lexical_cast<string>(*it);
          h2 = true;
        }
        else str_right+=","+lexical_cast<string>(*it);
    }
  of<<"\\fmfleft{"<<str_left<<"}\n";
  of<<"\\fmfright{"<<str_right<<"}\n";

  cout<<"work, left "<<str_left<< " right "<<str_right<<endl;
  print(G);
  /// Diagram output
  for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
    if((G[*ei].type == 1) && is_external_vert(vertex(G.global_to_local(G[*ei].to_vertex),G),G) && !G[*ei].visited)
      {
        print(G);
        G[G.global_to_local(vertex(G[*ei].to_vertex,G))].visited = true;
	G[*ei].visited = true;
	edis = *ei;
	bool visited_spinor_line = false;
        cout<<" work for,if"<<endl;
	while(!visited_spinor_line)
	  {
            //            of<<"\\fmf{fermion}{"<<G.global_to_local(G[edis].from_vertex) << "," << G.global_to_local(G[edis].to_vertex)<<"}\n";
            //    of<<"\\fmfv{l="<<vertex(G.global_to_local(G[edis].to_vertex),G)<<",l.a= 120,l.d=.02w}{"<<vertex(G[edis].to_vertex,G)<<"}\n";
            //of<<"\\fmfv{l="<<vertex(G.global_to_local(G[edis].from_vertex),G)<<",l.a= 120,l.d=.02w}{"<<vertex(G.global_to_local(G[edis].from_vertex),G)<<"}\n";

            // G[vertex(G.global_to_local(G[edis].from_vertex),G)].visited = true;
	    Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
            bool lng = false;
            //  for(tie(out_ei,out_ei_end) = out_edges(vertex(G.global_to_local(G[edis].from_vertex),G),G);out_ei != out_ei_end;++out_ei)
              if((G[*out_ei].type == 1) && (!G[*out_ei].visited)) 
                {
                  out_ei_tmp = out_ei;
                  lng = true;
                }
            if(lng == true)
              {
                G[*out_ei_tmp].visited = true;
                edis = *out_ei_tmp;
              }
            else visited_spinor_line = true;
	  }
      }
   // spinor edges in loops
  for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
    if((G[*ei].type == 1) && !G[*ei].visited)
      {
	bool visited_loop = false;
	edis = *ei;
	while(!visited_loop)
	  {
            //  G[vertex(G.global_to_local(G[edis].to_vertex),G)].visited = true;
        // of<<"\\fmf{fermion,left}{"<<G.global_to_local(G[edis].from_vertex) << "," << G.global_to_local(G[edis].to_vertex)<<"}\n";
            //   of<<"\\fmfv{l="<<vertex(G.global_to_local(G[edis].to_vertex),G)<<",l.a= 120,l.d=.02w}{"<<vertex(G.global_to_local(G[edis].to_vertex),G)<<"}\n";
            // of<<"\\fmfv{l="<<vertex(G.global_to_local(G[edis].from_vertex),G)<<",l.a= 120,l.d=.02w}{"<<vertex(G.global_to_local(G[edis].from_vertex),G)<<"}\n";
	    G[edis].visited = true;
            bool lng = false;
	    Graph::out_edge_iterator out_ei,out_ei_end,out_ei_tmp;
            // for(tie(out_ei,out_ei_end) = out_edges(vertex(G.global_to_local(G[edis].from_vertex),G),G); out_ei != out_ei_end; ++out_ei)
              if((G[*out_ei].type == 1) && (!G[*out_ei].visited))
                {
                  out_ei_tmp = out_ei;
                  lng = true;
                }
            if(lng == true)edis = *out_ei_tmp;		
            else visited_loop = true;
	      
	  }
      }
  of<<"\\fmffreeze\n";
  // other edges and other verticies
    for(tie(ei, ei_end) = edges(G); ei != ei_end; ++ei) 
    if(!G[*ei].visited)
      {
        // of<<"\\fmf{boson,left}{"<<G.global_to_local(G[*ei].from_vertex) << "," << G.global_to_local(G[*ei].to_vertex)<<"}\n";
        //of<<"\\fmfv{l="<<vertex(G.global_to_local(G[*ei].to_vertex),G)<<",l.a= 120,l.d=.02w}{"<<vertex(G.global_to_local(G[*ei].to_vertex),G)<<"}\n";
        //of<<"\\fmfv{l="<<vertex(G.global_to_local(G[*ei].from_vertex),G)<<",l.a= 120,l.d=.02w}{"<<vertex(G.global_to_local(G[*ei].from_vertex),G)<<"}\n";

	G[*ei].visited = true;
      }
  /*  for(tie(vi, vi_end) = vertices(G); vi != vi_end; ++vi)
    {
      if(G[*vi].visited)
        // G[*vi].visited = true;
        //out_degree(vdis,Gtmp)
      //      of<<"\\fmflabel{"<<*vi<<"}{"<<*vi<<"}\n";
        //of<<"\\fmfv{l="<<*vi<<",l.a= 120,l.d=.02w}{"<<*vi<<"}\n";
        }
  */
  of<<"\\end{fmfchar*}\n\\end{fmffile}\n";
  // cout<<"Drawing GRAPH no."<<get_property(G,graph_num)<<endl;
  return of.str();
  }
  else
    {
      string str("not a subgraph");
      return str;
    }
}



void afdcpp::fill_in_dia_ex()
{
  BOOST_FOREACH(Graph g1,gv11)
    {
      //print(g1);
      cout<<draw(g1)<<endl;
      in_dia_ex.append(ex_from_graph(g1));
    }
}

ex afdcpp::scalarize_vectors(ex in_ex)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | Convert indexed objects to scalars  | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
 scalar_products ind_sp;

  ind_sp.add(k1, p1, k1*p1);
  ind_sp.add(k1, p2, k1*p2);
  ind_sp.add(k1, q1, k1*q1);
  ind_sp.add(k1, q2, k1*q2);
  ind_sp.add(k2, p1, k2*p1);
  ind_sp.add(k2, p2, k2*p2);
  ind_sp.add(k2, q1, k2*q1);
  ind_sp.add(k2, q2, k2*q2);
  ind_sp.add(k3, p1, k3*p1);
  ind_sp.add(k3, p2, k3*p2);
  ind_sp.add(k3, q1, k3*q1);
  ind_sp.add(k3, q2, k3*q2);
  ind_sp.add(k4, p1, k4*p1);
  ind_sp.add(k4, p2, k4*p2);
  ind_sp.add(k4, q1, k4*q1);
  ind_sp.add(k4, q2, k4*q2);
  ind_sp.add(k5, p1, k5*p1);
  ind_sp.add(k5, p2, k5*p2);
  ind_sp.add(k5, q1, k5*q1);
  ind_sp.add(k5, q2, k5*q2);
  ind_sp.add(k6, p1, k6*p1);
  ind_sp.add(k6, p2, k6*p2);
  ind_sp.add(k6, q1, k6*q1);
  ind_sp.add(k6, q2, k6*q2);



  ind_sp.add(k1, k1, k1*k1);
  ind_sp.add(k1, k2, k1*k2);
  ind_sp.add(k1, k3, k1*k3);
  ind_sp.add(k1, k4, k1*k4);
  ind_sp.add(k1, k5, k1*k5);
  ind_sp.add(k1, k6, k1*k6);
  ind_sp.add(k2, k2, k2*k2);
  ind_sp.add(k2, k3, k2*k3);
  ind_sp.add(k2, k4, k2*k4);
  ind_sp.add(k2, k5, k2*k5);
  ind_sp.add(k2, k6, k2*k6);
  ind_sp.add(k3, k3, k3*k3);
  ind_sp.add(k3, k4, k3*k4);
  ind_sp.add(k3, k5, k3*k5);
  ind_sp.add(k3, k6, k3*k6);
  ind_sp.add(k4, k4, k4*k4);
  ind_sp.add(k4, k5, k4*k5);
  ind_sp.add(k4, k6, k4*k6);
  ind_sp.add(k5, k5, k5*k5);
  ind_sp.add(k5, k6, k5*k6);
  ind_sp.add(k6, k6, k6*k6);

  ind_sp.add(p1,p1, p1*p1);
  ind_sp.add(p2,p2, p2*p2);
  ind_sp.add(q1,q1, q1*q1);
  ind_sp.add(q2,q2, q2*q2);
  ind_sp.add(p1,q1, p1*q1);
  ind_sp.add(p2,q2, p1*q2);
  ind_sp.add(p2,q1, p2*q1);
  ind_sp.add(p2,q2, p2*q2);
  ind_sp.add(p1,p2, p1*p2);
  ind_sp.add(q1,q2, q1*q2);
  return in_ex.expand(expand_options::expand_indexed).simplify_indexed(ind_sp);
}

ex afdcpp::scalarize_l_vectors(ex in_ex)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | NEW procedure for L vectors         | ///
/// | Convert indexed objects to scalars  | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
 scalar_products ind_sp;

  ind_sp.add(k1, l1, k1*l1);
  ind_sp.add(k1, l2, k1*l2);
  ind_sp.add(k1, l3, k1*l3);
  ind_sp.add(k1, l4, k1*l4);
  ind_sp.add(k2, l1, k2*l1);
  ind_sp.add(k2, l2, k2*l2);
  ind_sp.add(k2, l3, k2*l3);
  ind_sp.add(k2, l4, k2*l4);
  ind_sp.add(k3, l1, k3*l1);
  ind_sp.add(k3, l2, k3*l2);
  ind_sp.add(k3, l3, k3*l3);
  ind_sp.add(k3, l4, k3*l4);
  ind_sp.add(k4, l1, k4*l1);
  ind_sp.add(k4, l2, k4*l2);
  ind_sp.add(k4, l3, k4*l3);
  ind_sp.add(k4, l4, k4*l4);
  ind_sp.add(k5, l1, k5*l1);
  ind_sp.add(k5, l2, k5*l2);
  ind_sp.add(k5, l3, k5*l3);
  ind_sp.add(k5, l4, k5*l4);
  ind_sp.add(k6, l1, k6*l1);
  ind_sp.add(k6, l2, k6*l2);
  ind_sp.add(k6, l3, k6*l3);
  ind_sp.add(k6, l4, k6*l4);



  ind_sp.add(k1, k1, k1*k1);
  ind_sp.add(k1, k2, k1*k2);
  ind_sp.add(k1, k3, k1*k3);
  ind_sp.add(k1, k4, k1*k4);
  ind_sp.add(k1, k5, k1*k5);
  ind_sp.add(k1, k6, k1*k6);
  ind_sp.add(k2, k2, k2*k2);
  ind_sp.add(k2, k3, k2*k3);
  ind_sp.add(k2, k4, k2*k4);
  ind_sp.add(k2, k5, k2*k5);
  ind_sp.add(k2, k6, k2*k6);
  ind_sp.add(k3, k3, k3*k3);
  ind_sp.add(k3, k4, k3*k4);
  ind_sp.add(k3, k5, k3*k5);
  ind_sp.add(k3, k6, k3*k6);
  ind_sp.add(k4, k4, k4*k4);
  ind_sp.add(k4, k5, k4*k5);
  ind_sp.add(k4, k6, k4*k6);
  ind_sp.add(k5, k5, k5*k5);
  ind_sp.add(k5, k6, k5*k6);
  ind_sp.add(k6, k6, k6*k6);

  ind_sp.add(l1,l1, l1*l1);
  ind_sp.add(l2,l2, l2*l2);
  ind_sp.add(l3,l3, l3*l3);
  ind_sp.add(l4,l4, l4*l4);
  ind_sp.add(l1,l3, l1*l3);
  ind_sp.add(l2,l4, l1*l4);
  ind_sp.add(l2,l3, l2*l3);
  ind_sp.add(l2,l4, l2*l4);
  ind_sp.add(l1,l2, l1*l2);
  ind_sp.add(l3,l4, l3*l4);
  return in_ex.expand(expand_options::expand_indexed).simplify_indexed(ind_sp);
}

///please mark all momentums before input,otherwise new moments are input only!!!
ex afdcpp::subs_ext_mom(ex in_ex,exvector ext_vec)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | External momentum substitution      | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  //symbol l1("l1"),l2("l2"),l3("l3"),l4("l4"),l5("l5");
  //const ex lm[] = {l1,l2,l3,l4,l5};
  int n = 0;
  //std::list<ex> m_lst_lst;
  //list of all momentum symbols in external lines 
  std::list<symbol> m_lst;
  BOOST_FOREACH(ex e478,ext_vec)
    {

      cout<<"mom of vert: "<<endl;
      std::copy(e478.preorder_begin(),e478.preorder_end(),std::ostream_iterator<ex>(cout, "\n"));
      cout<<"my way:"<<endl;
      
      for (const_preorder_iterator it = e478.preorder_begin();it!=e478.preorder_end();++it)
        {
          if(is_a<symbol>(*it))m_lst.push_back(ex_to<symbol>(*it));
          cout<<"---> "<<(*it)<<" is a symbol: "<<is_a<symbol>(*it)<<endl;
        }
        
          //      cout<<"solution subs "<<(*it)<<" = "<<lsolve(e478==lm[n],(*it))<<endl;
    }
  BOOST_FOREACH(symbol e634,m_lst)
    {
      cout<<"not unique "<<e634<<endl;
    }
  m_lst.unique();
  BOOST_FOREACH(symbol e634,m_lst)
    {
      cout<<"unique "<<e634<<endl;
    }
  /// creation of momentum index map
  //std::map<symbol,int> mom_idx;
  exhashmap<int> mom_idx;
  BOOST_FOREACH(ex e789,m_lst)
    {
      BOOST_FOREACH(ex e45,ext_vec)
        {
          if(e45.has(e789))
            {
              cout<<e45<<" has "<<e789<<endl;
              mom_idx[ex_to<symbol>(e789)]+=1;
            }
        }
    }
  for(exhashmap<int>::iterator it = mom_idx.begin();it != mom_idx.end();++it)
    cout<<"mom: "<<it->first<<" mom_idx: "<<it->second<<endl;
  /// finding of minimal momentum index in each ext_mom
  /// and mom map creation
  //list in form (ext_mom,min),where min is a momentum with minimal index
  std::vector< std::pair<ex,symbol> > em_m_lst;
  BOOST_FOREACH(ex e45,ext_vec)
    {
      symbol buf_mom;
      int buf_idx = 1000;
      BOOST_FOREACH(symbol e789,m_lst)
        {
          if(e45.has(e789) && (mom_idx[e789] < buf_idx))
            {
              buf_idx = mom_idx[e789];
              buf_mom = e789;
            }
        }
      em_m_lst.push_back(make_pair(e45,buf_mom));
    }
  for(std::vector<std::pair<ex,symbol> >::iterator i = em_m_lst.begin();i!=(em_m_lst.end()-1);++i)
    if(mom_idx[i->second]>mom_idx[(i+1)->second])swap(*i,*(i+1));
  typedef std::pair<ex,symbol> mompair;
  
  exmap mom_subs_map;    
  BOOST_FOREACH(mompair p1,em_m_lst)
    {
      mom_subs_map[p1.second] = lsolve(p1.first==lm[n],p1.second);
      cout<<"Ext mom: "<<p1.first<<" min "<<p1.second<< " index "<<mom_idx[p1.second]<<endl;
      n++;
    }
  cout<<colors::blue<<"SMOM:  "<<in_ex.subs(mom_subs_map)<<colors::def_col<<endl;
  return in_ex.subs(mom_subs_map);
}

void afdcpp::feynman(ex inex,exvector loop_vec,exvector ext_vec)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | Integral Feynman parametrization    | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  const int n_loops = 3;
  /// list of loop moments
  lst k_lst;
  k_lst = k1,k2,k3,k4,k5,k6;
  /// list of loop moment scalar products
  lst kk_lst;
  kk_lst = k1_k1,
    k1_k2,
    k1_k3,
    k1_k4,
    k1_k5,
    k1_k6,
    k2_k2,
    k2_k3,
    k2_k4,
    k2_k5,
    k2_k6,
    k3_k3,
    k3_k4,
    k3_k5,
    k3_k6,
    k4_k4,
    k4_k5,
    k4_k6,
    k5_k5,
    k5_k6,
    k6_k6;

  lst l;
  l = k1_p1,
    k1_p2,
    k1_q1,
    k1_q2,
    k2_p1,
    k2_p2,
    k2_q1,
    k2_q2,
    k3_p1,
    k3_p2,
    k3_q1,
    k3_q2,
    k4_p1,
    k4_p2,
    k4_q1,
    k4_q2,
    k5_p1,
    k5_p2,
    k5_q1,
    k5_q2,
    k6_p1,
    k6_p2,
    k6_q1,
    k6_q2,
    k1_k1,
    k1_k2,
    k1_k3,
    k1_k4,
    k1_k5,
    k1_k6,
    k2_k2,
    k2_k3,
    k2_k4,
    k2_k5,
    k2_k6,
    k3_k3,
    k3_k4,
    k3_k5,
    k3_k6,
    k4_k4,
    k4_k5,
    k4_k6,
    k5_k5,
    k5_k6,
    k6_k6;


  
  // inex = inex.subs(mS == 0); 
  //  inex = inex.expand(expand_options::expand_indexed).simplify_indexed(kk_lst);
  /// Preparing NUMERATOR
  ///////////////////////////////////
  /// Convert to Rational expression/
  inex = factor(inex).normal();//////
  ///////////////////////////////////
  ex N = inex.numer();
  cout<<"NUMER wo TRACE: "<< N<<endl<<endl;
  //ex sN = N.expand(expand_options::expand_indexed).simplify_indexed(kk_lst);
  //ex sN = sN.collect(l);
  cout<<endl;
  cout<<"NUMERATOR "<<factor(N)<<"expr"<<inex<<endl;
  /// NUMERATOR prepared
  
  /// Preparing DENOMINATOR
      /// DENUMERATOR  propogators for solving
  lst xp; //coefficients in propagators list
  ex D = inex.denom();
  //      std::cout<<"D: "<<inex.denom()<<std::endl;
  //      D.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
  ex feyn_d;
  unsigned int n = 0;
  string  str;
  ex  sD;   // Linear combination of propagators
  lst prop; // list of propagators
  lst Pr;   // list of P_i symbols
  ex  Px;   // Denominator in form P_i*X_i
  // symbol s("x_0"); // addition of free term
  // xp.append(s);    //
  //sD = s;          //
  //Px = s;          // P_0==1 (Scalar term)
  for(const_iterator i = D.begin(); i != D.end(); ++i)
    {
      if( !is_a<GiNaC::numeric>(*i))
        {
          cout<<"subs: "<<factor( (*i).expand(expand_options::expand_indexed).simplify_indexed(global_sp())) <<endl;
          
          cout<<"Dividing: "<<N.collect(*i)<<endl;
          str = "x_" + lexical_cast<string>(n);
          symbol s(str);
          xp.append(s);
          sD += xp[n] * (*i); 
          prop.append(*i);
          str = "P" + lexical_cast<string>(n);
          symbol p_sym(str);
          Pr.append(p_sym);
          Px += s*p_sym;
        }
      ++n;
    }
  
  cout<<"Number of props: "<<xp.nops()<<" P_i*X_i = "<<Px<<endl;
  /// from indexed to scalar expressions:
  ///
  /// sD prepared for matrix building
  sD = scalarize_l_vectors(sD);
  //sD = sD.expand(expand_options::expand_indexed).simplify_indexed(ind_sp);
  cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
  cout<<sD<<endl<<"^^^^^^^^^^^^^^^^^^^^^^^^^^^"<<endl;
  if(sD.has(k1))cout<<"HAS k1"<<endl;
  cout<<"Coeff k1:"<<sD.coeff(k2*k3,1)<<endl;
  /// DENOMINATOR Prepared
  ex sDtmp = sD;// expr for modification
  matrix M(loop_vec.size(),loop_vec.size());
  matrix Q(loop_vec.size(),1);//column
  GiNaC::numeric half(1,2);
  for(exvector::iterator itr = loop_vec.begin();itr!=loop_vec.end();++itr)
    for(exvector::iterator itc = itr;itc!=loop_vec.end();++itc)
      if(itr == itc)
        {
          M(int(itr - loop_vec.begin()),int(itc - loop_vec.begin())) = -1*sD.coeff((*itr)*(*itc),1);
          cout<<(*itr)*(*itc)<<"M("<<int(itr - loop_vec.begin())<<","<<int(itc - loop_vec.begin())<<") coeff "<<sD.coeff((*itr)*(*itc),1)<<endl;
          sDtmp -= (*itr)*(*itc)*sD.coeff((*itr)*(*itc),1);
        }
      else
        {
          M(int(itr - loop_vec.begin()),int(itc - loop_vec.begin())) = -1*half*sD.coeff((*itr),1).coeff((*itc),1);
          M(int(itc - loop_vec.begin()),int(itr - loop_vec.begin())) = -1*half*sD.coeff((*itr),1).coeff((*itc),1);
          cout<<(*itr)*(*itc)<<"M("<<int(itr - loop_vec.begin())<<","<<int(itc - loop_vec.begin())<<") coeff "<<sD.coeff((*itr),1).coeff((*itc),1)<<endl;
          sDtmp -= (*itr)*(*itc)*sD.coeff((*itr),1).coeff((*itc),1);
        }
  cout<<colors::yelow<<"M: "<<M<<colors::def_col<<endl;
  sDtmp = sDtmp.expand();
  cout<<"Expr linear on external momentum: "<<sDtmp.expand()<<endl;

  for(exvector::iterator itr = loop_vec.begin();itr!=loop_vec.end();++itr)
    {
      Q(int(itr - loop_vec.begin()),0) = half*sDtmp.coeff((*itr),1);
      sDtmp -= (*itr)*sDtmp.coeff((*itr),1);
    }
  cout<<colors::yelow<<"Q: "<<Q<<colors::def_col<<endl;
  sDtmp = sDtmp.expand();
  ex minusJ = sDtmp;
  cout<<colors::yelow<<"-J: "<<minusJ<<colors::def_col<<endl;
  ex U = M.determinant();
  ex F = expand(M.determinant()*(minusJ+Q.transpose().mul(M.inverse()).mul(Q)(0,0)));
  cout<<"U= "<<U<<endl<<"F= "<<F<<endl;
  //Zero momentum
  exmap zerosubs;
  BOOST_FOREACH(ex e478,ext_vec)
    {
      cout<<"mom of vert: "<<endl;
      std::copy(e478.preorder_begin(),e478.preorder_end(),std::ostream_iterator<ex>(cout, "\n"));
      cout<<"my way:"<<endl;
      for(const_preorder_iterator it = e478.preorder_begin();it!=e478.preorder_end();++it)
        {
          cout<<"---> "<<(*it)<<" is a symbol: "<<is_a<symbol>(*it)<<endl;
          if(is_a<symbol>(*it))zerosubs[(*it)] = 0.000001;
        }
      cout<<zerosubs<<" end of my way"<<endl;

    }
  cout<<"U= "<<U.subs(zerosubs)<<endl<<"F= "<<F.subs(zerosubs)<<endl;
  U = U.subs(lst(l1==-1,l2==1,l3==1,l4==1));
  F = F.subs(lst(l1==-1,l2==1,l3==1,l4==1)).subs(mS == 1);
  try{
    using namespace sector_decomposition;
    exvector poly_list(2);
    poly_list[0] = U;
    poly_list[1] = F;
    std::vector<exponent> c(2);
    cout<<"exp test "<<(loop_vec.size()+1)<<endl;
    c[0] = exponent(xp.nops()-2*(loop_vec.size()+1),loop_vec.size()+1);
    c[1] = exponent(-xp.nops()+2*loop_vec.size(),-loop_vec.size());
     /// construct vector of feynman parameters
    exvector parameters;
    for(lst::const_iterator i = xp.begin();i!=xp.end();i++)
      parameters.push_back((*i));

  cout<<parameters<<endl;
    symbol eps("eps");
    integrand my_integrand(xp.nops(),poly_list,c);
    monte_carlo_parameters mc_p(5,15,10000,100000);
    for(int order=-2*loop_vec.size();order <=0;order++)
      {
      integration_data gdata(parameters,eps,order);
      monte_carlo_result res = do_sector_decomposition(gdata,my_integrand,mc_p,5);
      std::cout << "Order " << pow(eps,order) << ": " << res.get_mean() 
                << " +/- " << res.get_error() << std::endl;
      
      }
  } catch (std::exception &e)
    {
      std::cout << "Exception : " << e.what() << std::endl;
    }
  

}

/// PROCEDURE FOR ONELOOP NUMERATOR EXPR IN TERMS OF PROPS
void afdcpp::reduze_out(ex inex)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  try
    {
      lst l;
      l = k1_p1,
        k1_p2,
        k1_q1,
        k1_q2,
        k2_p1,
        k2_p2,
        k2_q1,
        k2_q2,
        k3_p1,
        k3_p2,
        k3_q1,
        k3_q2,
        k4_p1,
        k4_p2,
        k4_q1,
        k4_q2,
        k5_p1,
        k5_p2,
        k5_q1,
        k5_q2,
        k6_p1,
        k6_p2,
        k6_q1,
        k6_q2,
        k1_k1,
        k1_k2,
        k1_k3,
        k1_k4,
        k1_k5,
        k1_k6,
        k2_k2,
        k2_k3,
        k2_k4,
        k2_k5,
        k2_k6,
        k3_k3,
        k3_k4,
        k3_k5,
        k3_k6,
        k4_k4,
        k4_k5,
        k4_k6,
        k5_k5,
        k5_k6,
        k6_k6;
      
      inex = inex.subs(q1 == p1); 
      inex = inex.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
      /// Preparing NUMERATOR
      ///////////////////////////////////
      /// Convert to Rational expression/
      inex = factor(inex).normal();//////
      ///////////////////////////////////
      ex N = inex.numer();
      cout<<"NUMER wo TRACE: "<< N<<endl<<endl;
      ex sN = dirac_trace(N).expand(expand_options::expand_indexed).simplify_indexed(global_sp());
      //      sN = sN.collect(l);
      cout<<endl;
      cout<<"NUMERATOR "<<factor(sN)<<endl;
      /// NUMERATOR prepared

      /// Preparing DENOMINATOR
      /// DENUMERATOR  propogators for solving
      lst xp; //coefficients in propagators list
      ex D = inex.denom();
      //      std::cout<<"D: "<<inex.denom()<<std::endl;
      //      D.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
      ex feyn_d;
      unsigned int n = 1;
      string  str;
      ex  sD;   // Linear combination of propagators
      lst prop; // list of propagators
      lst Pr;   // list of P_i symbols
      ex  Px;   // Denominator in form P_i*X_i
      symbol s("x_0"); // addition of free term
      xp.append(s);    //
      sD = s;          //
      Px = s;          // P_0==1 (Scalar term)
      for(const_iterator i = D.begin(); i != D.end(); ++i)
        {
          if( !is_a<GiNaC::numeric>(*i))
            {
              cout<<"subs: "<<factor( (*i).expand(expand_options::expand_indexed).simplify_indexed(global_sp())) <<endl;
              
              cout<<"Dividing: "<<N.collect(*i)<<endl;
              str = "x_" + lexical_cast<string>(n);
              symbol s(str);
              xp.append(s);
              sD += xp[n] * (*i); 
              prop.append(*i);
              str = "P" + lexical_cast<string>(n);
              symbol p_sym(str);
              Pr.append(p_sym);
              Px += s*p_sym;
            }
          ++n;
        }
          
      cout<<" P_i*X_i = "<<Px<<endl;
      /// sD prepared for matrix building
      sD = sD.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
      cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
      cout<<sD<<endl<<"^^^^^^^^^^^^^^^^^^^^^^^^^^^"<<endl;
      /// DENOMINATOR Prepared
      
      ex N_in_prop;
      lst sp_in_N;
      exmap sp_prop_map;


      for (lst::const_iterator it = l.begin(); it != l.end(); ++it)
        if(sN.has(*it))sp_in_N.append(*it);

      cout<<"N has"<<sp_in_N<<"dfdf"<<endl;

      for(lst::const_iterator itsp = sp_in_N.begin(); itsp != sp_in_N.end(); ++itsp)
        {
          ex sDtmp;          
          lst eqns;
          cout<<"TERM: "<<(*itsp)<<endl<<endl;
          for (lst::const_iterator it1 = l.begin(); it1 != l.end(); ++it1)
            {
              if(sD.has(*it1))
                {
                  if((*itsp) ==  *it1)
                    {
                      eqns.append(sD.coeff(*it1,1) == 1);
                      sDtmp += sD.coeff(*it1,1)*(*it1);
                      cout<<"\n Founded dividing: "<<sD.coeff(*it1,1)*(*it1)<<endl;
                    }
                  else
                    {
                      eqns.append(sD.coeff(*it1,1) == 0);
                      sDtmp += sD.coeff(*it1,1)*(*it1);
                      cout<<"\n NOT  dividing: "<<sD.coeff(*it1,1)*(*it1)<<endl;
                    }

                 
                  //cout<<endl<<"SD value:"<<sD<<endl;
                  //     cout<<"SD value end ^^"<<endl<<endl;
                }
               
            }
          //          cout<<"SD "<<sD<<endl<<endl;
          cout<<"SD - SDtmp"<<(sD-sDtmp).collect(l)<<endl;
          ex sDx0 = (sD-sDtmp).collect(l);
          eqns.append(0==sDx0);
          
          cout<<eqns.nops()<<" equations with "<<xp.nops() <<" variables to be solved: " <<eqns<<endl;
          cout<<"solved "<<lsolve(eqns,xp)<<endl;
          lst null_terms;
                                        
          for (lst::const_iterator it1 = xp.begin(); it1 != xp.end(); ++it1)
            {
              null_terms.append((*it1)==0);
            }

          sp_prop_map[(*itsp)] = Px.subs(lsolve(eqns,xp)).subs(null_terms);
          //          sp_in_prop.append(Px.subs(lsolve(eqns,xp))+xp[0];
        }

      std::cout<<std::endl<<"##############################"<<std::endl;
      std::cout<<"## eqs: "<<sp_prop_map<<std::endl;
      std::cout<<std::endl<<"##############################"<<std::endl;

      ex Ntmp;
      Ntmp = expand(sN.subs(sp_prop_map,subs_options::algebraic));
      /// NUMERATOR IN TERMS of P_i
      ex Ncol;
      Ncol = factor(Ntmp).expand().collect(Pr);
      //      Ncol = Ntmp.collect(Pr);
      cout<<"\nNUMERATOR in terms of P_i: "<<Ncol<<endl;
      cout<<"\nNumber of terms in N "<<Ncol.nops()<<endl;


    }
  catch(std::exception &p)
    {
      std::cout<<std::endl<<"##############################"<<std::endl;
      std::cout<<"## ERROR: "<<p.what()<<std::endl;
      std::cout<<std::endl<<"##############################"<<std::endl;
    }
}

void afdcpp::reduze_in_dia_ex()
{
  //  BOOST_FOREACH(ex e1,in_dia_ex)
    {
      reduze_out(in_dia_ex.op(1));

      divrat(in_dia_ex.op(1),p1);
    }
}

void afdcpp::qed_renorm()
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |Example of qed diaram renormalization| ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  //Graph_1PI_lst tmplst = make_subgraphs(gv11[0]);
  //get_terms(tmplst);
}


void afdcpp::scalar_renorm()
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |Example of qed diaram renormalization| ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  // Graph_1PI_lst tmplst = make_subgraphs(gv11[0]);
  //ex substr = get_terms(tmplst);
  //ex G_ex = ex_from_graph(gv11[0]);
 // feynman(G_ex);
}

void afdcpp::init()
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// | Init after constructor execution    | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
try
  {
  //gsp = global_sp();
  //fill_in_dia_ex();
  // reduze_in_dia_ex();
  //M(gv11[0]);

  ///disable cerr!!!!!

  // std::stringstream output;
  //redirect cerr output
  //freopen ("errfile.txt","w",stderr);


  cout<<"inited\n";
  cout <<std::setw (10);
  cout << 77 <<std::setw (10)<<88<<endl ;
  print_graph(gv11[0]);
  logo();
  p_table();
  ex in_ex_gr = ex_from_graph(gv11[0]);
  //cout<<"Tr[N] "<<endl<<dirac_trace(in_ex_gr.numer()).expand(expand_options::expand_indexed).simplify_indexed(global_sp());
  //feynman(in_ex_gr);
  //  qed_renorm();
    
  //scalar_renorm();
///works part[
  //raph_1PI_lst tmplst = make_subgraphs(gv11[0]);
  //et_terms(tmplst);
  ///]


//draw_lst(tmplst);
  print(gv11[0]);
  Graph g67(gv11[0]);
  //!!!!!  renorm_dia d23(g67,ext_legs);
  //  cout<<draw2(gv11[0]);
  int n = 1;


  BOOST_FOREACH(Graph g,gv11)
    {
       renorm_dia d23(g,ext_legs);
      //cout<<n<<endl;
      //print_graph(g);
      //if(n == gv11.size())draw2(g,true);
      //else draw2(g,false);
      dr1.add(g);
      //cout<<"-----------------------------------------------"<<endl;
           n++;
    }
  //dr1.print();
  
#ifndef BOOST_GRAPH_NO_BUNDLED_PROPERTIES
  cout<<"exists"<<endl;
#endif
  }catch(std::exception &p)
    {
      cout<<"******************************************************"<<endl;
      cout<<">>>  ERROR in class \"AFDCPP::init\" :\n |___>   "<<setw(40)<<std::internal<<p.what()<<endl;
      cout<<"******************************************************"<<endl;
      //      throw std::logic_error(std::string("In function AFDCPP::INIT:\n |___> ")+p.what());
     }

}


/// END of Constructing program Functions 
ex afdcpp::EeA(varidx& index)
{
  ex e=dirac_gamma(index);
  return e;
}

ex afdcpp::BbZ(varidx& index)
{
  ex e=dirac_gamma(index);
  e = e.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
  cout<<"\nDEBUG BbZ expr = " <<e.simplify_indexed()<<endl<<endl;

  return e;
}

ex afdcpp::BbGZ()
{
  ex e=lB;
  return e;
}

ex afdcpp::roZ(varidx& index, ex& mom1, ex& mom2, ex& mom3)
{
  try
    {
      ex e=indexed(mom1,index);
      e = e.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
      cout<<"\nDEBUG roZ expr = " <<e.simplify_indexed()<<endl<<endl;

      return e;
    }
  catch(std::exception &p)
    {
      std::cout<<std::endl<<"##############################"<<std::endl;
      std::cout<<"## ERROR: "<<p.what()<<std::endl;
      std::cout<<std::endl<<"##############################"<<std::endl;
    }
}

///Spinor propagators
ex afdcpp::S(symbol& sim,ex& mom)
{
  varidx tmp_idx1(symbol("dnu1"),4);
  //ex e = (dirac_slash(mom,4) + sim*dirac_ONE())/(pow(mom,2) - pow(sim,2));
  ex e = (dirac_slash(mom,4) + sim*dirac_ONE())/(indexed(mom,tmp_idx1)*indexed(mom,tmp_idx1.toggle_variance()) - pow(sim,2));
  //  ex e = (dirac_slash(mom,4) + sim*dirac_ONE())/((dirac_slash(mom,4) - sim*dirac_ONE())*(dirac_slash(mom,4) + sim*dirac_ONE()));
  return e;
}

ex afdcpp::B(symbol& sim,ex& mom)
{
  try
    {
      varidx tmp_idx1(symbol("dnu1"),4);
      ex e = (dirac_slash(mom,4) + sim*dirac_ONE())/(indexed(mom,tmp_idx1)*indexed(mom,tmp_idx1.toggle_variance()) - pow(sim,2));

      e = e;//.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
      cout<<"\nDEBUG B expr = " <<e.simplify_indexed()<<endl<<endl;

      return e;
    }
  catch(std::exception &p)
    {
      std::cout<<std::endl<<"##############################"<<std::endl;
      std::cout<<"## ERROR: "<<p.what()<<std::endl;
      std::cout<<std::endl<<"##############################"<<std::endl;
    }
}

/// Gauge bosons propagators 
/// All indices in propagators are covariant
ex afdcpp::D(symbol& sim,ex& mom,varidx& v_to,varidx& v_from)
{
  varidx tmp_idx1(symbol("dnu1"),4);
  ex e = lorentz_g(v_to.toggle_variance(),v_from.toggle_variance())/
    (indexed(mom,tmp_idx1)*indexed(mom,tmp_idx1.toggle_variance())-pow(sim,2));
  return e;
}

ex afdcpp::Z(symbol& sim,ex& mom,varidx& v_to,varidx& v_from)
{
  try
    {
      varidx tmp_idx1(symbol("dnu2"),4);
      varidx tmp_idx2(symbol("dnu3"),4);
      ex e3 =indexed(mom,tmp_idx2)*indexed(mom,tmp_idx2.toggle_variance());
      //    e3 = e3.expand(expand_options::expand_indexed);
      cout<<e3.simplify_indexed(global_sp())<<endl;
    
        
      ex e0 = ( lorentz_g(v_to.toggle_variance(),v_from.toggle_variance())*(indexed(mom,tmp_idx2)*indexed(mom,tmp_idx2.toggle_variance())-xiZ*pow(sim,2)) -
                indexed(mom,v_to.toggle_variance()) * indexed(mom,v_from.toggle_variance())*(1-xiZ));
      //e0 = e0.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
      ex e1 = (indexed(mom,tmp_idx2)*indexed(mom,tmp_idx2.toggle_variance())-xiZ*pow(sim,2));
            
      //e1 = e1.expand(expand_options::expand_indexed).simplify_indexed(global_sp());            
    
      ex e2 = (indexed(mom,tmp_idx1)*indexed(mom,tmp_idx1.toggle_variance())-pow(sim,2));
      //e2 = e2.expand(expand_options::expand_indexed).simplify_indexed(global_sp());

        
    
      ex e = e0/e1/e2;
      //      e = e.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
      cout<<"\nDEBUG Z expr = " <<e.simplify_indexed()<<endl<<endl;
      return e;
    }
  catch(std::exception &p)
    {
      std::cout<<std::endl<<"##############################"<<std::endl;
      std::cout<<"## ERROR: "<<p.what()<<std::endl;
      std::cout<<std::endl<<"##############################"<<std::endl;
    }
}

/// Goldstone bosons propagators
ex afdcpp::GZ(symbol& sim,ex& mom)
{
  try{
    varidx tmp_idx8(symbol("dnu8"),4);
    ex e = 1/(indexed(mom,tmp_idx8)*indexed(mom,tmp_idx8.toggle_variance()) - pow(sim,2));//.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
    
    //          e = e.expand(expand_options::expand_indexed).simplify_indexed(global_sp());
    cout<<"\nDEBUG GZ expr = " <<e.simplify_indexed()<<endl<<endl;

    return e;
  }
  catch(std::exception &p)
    {
      std::cout<<std::endl<<"##############################"<<std::endl;
      std::cout<<"## ERROR: "<<p.what()<<std::endl;
      std::cout<<std::endl<<"##############################"<<std::endl;
    }
}

/// Prepare diagram expression for calculation divergeces
/// Trace calculation, external momentums substituting, etc.


/// R-operation part of code

/// Differentiation propagator type expression
ex afdcpp::divpoly(ex ine,ex dp,varidx d_idx)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
try
  {
    if (ine.has(dp)) /// div if expr depends on dp ONLY!!!
      {
        if(is_a<power>(ine))
          {
            
            power p_diff(ine.op(0),ine.op(1)-1);
            return ine.op(1)*p_diff*divpoly(ine.op(0),dp,d_idx);
          }
        else if(is_a<mul>(ine)) /// If it's a mul at first look
        {
          ine.expand();
          ex diff_out,mult_wo_diff;
          diff_out = 0;
         
          cout<<"\n\nDiff Multyplication function...\n "<<ine<<endl;
          for(const_iterator i = ine.begin(); i != ine.end(); ++i)
            {
              if((*i).has(dp))
                {
                  //cout<<std::setw(40)<<(*i)<<endl;
                  mult_wo_diff = 1;
                  for(const_iterator j = ine.begin(); j != ine.end(); ++j)
                    if (j != i) mult_wo_diff *= (*j);
                  //cout<< "mul constructed"<<endl;
                  //ex tmp_ex = (*i);
                  //tmp_ex = divpoly(tmp_ex);
                  diff_out += (mult_wo_diff *divpoly(*i,dp,d_idx));
                }
            }
          cout<<"Before mul diff "<<ine<<"\nAfter diff mul func "<<diff_out<<endl;
          return diff_out;
        }

      else if(is_a<add>(ine))  /// If it's a sum at first look
        {
          ex diff_out = 0;
          cout<<"\n\nDiff Addition function...\n "<<ine<<endl;
          for(const_iterator i = ine.begin(); i != ine.end(); ++i)
            if((*i).has(dp))
              diff_out += divpoly(*i,dp,d_idx);
          cout<<"Before add diff"<<ine<<"\nAfter add diff "<<diff_out<<endl;
          return diff_out;
        }


      /////////////////////////////////////////////////////////////
      /////                IMPORTANT!!!                     ///////
      /////differentiating over kontravariant index!!! p~mu ///////
      /////////////////////////////////////////////////////////////
      else if(is_a<indexed>(ine)) /// assuming that indexed may be only the sum
        {
          cout<< "Diff indexed "<<ine;
          exmap idx_fnd;
          if(ine.match(indexed(dp,varidx(wild(),4,true)),idx_fnd))
            {
              cout<<"Covariant simple expr "<<ine<<" with idx "<<idx_fnd[wild()]<<endl;
              /// d(p.mu)/d(p~nu) = g.mu.nu 
              exvector iv = ine.get_free_indices();
              if(iv.size() > 1) throw 13; /// Exception, more then one index in object
              cout<<" and diff is: "<<lorentz_g(iv[0],d_idx.toggle_variance())<<endl;              
              return lorentz_g(iv[0],d_idx.toggle_variance());                            
            }
          else if(ine.match(indexed(dp,varidx(wild(),4,false)),idx_fnd))
            {
              cout<<"Contravariant simple expr "<<ine<<" with idx "<<idx_fnd[wild()]<<endl;
              
              /// d(p~mu)/d(p~nu) = delta~mu.nu
              exvector iv = ine.get_free_indices(); 
              if(iv.size() > 1) throw 13; /// Exception, more then one index in object
              const varidx & vi = ex_to<varidx>(iv[0]);
              if (vi.is_contravariant()) 
                {
                  cout<<" and diff is: "<<delta_tensor(iv[0],d_idx.toggle_variance())<<endl;
               return delta_tensor(iv[0],d_idx.toggle_variance());
                }
              if (vi.is_covariant())
                {
                  cout<<" and diff is: "<<delta_tensor(iv[0],d_idx)<<endl;
               return delta_tensor(iv[0],d_idx);
                }
              //return delta_tensor(vi,d_idx.toggle_variance());
            }
          else
            {
              cout<<"\nFound default: "<<ine<<endl<<endl;
              cout<<"default and diff is: "<<divpoly(ine.expand(expand_options::expand_indexed),dp,d_idx)<<endl;
              return divpoly(ine.expand(expand_options::expand_indexed),dp,d_idx);
            }
        }
      else
        {
          cout<<std::setw(20)<<"Unknown expr"<<setw(20)<<ine.expand(expand_options::expand_indexed)<<endl;
          
          return ine;
        }
 
      }
  }
 catch(int n)
   {
     if(n ==13)
       {
         std::cout<<std::endl<<"##############################"<<std::endl;
         std::cout<<"## ERROR: "<<" More then one index in indexed object"<<std::endl;
         std::cout<<std::endl<<"##############################"<<std::endl;
       }
   }
}

void afdcpp::divrat(ex ine,ex dp)
{
  varidx to_dif(symbol("d_idx"),4);
  cout <<"DIV STARTED"<<endl;
  ine = ine.subs(q1 == p1);

for(const_iterator i = ine.begin(); i != ine.end(); ++i)
  {
    cout<<"term:  "<<(*i)<<" is a Power "<<is_a<power>(*i)<<endl;
    if(is_a<power>(*i)) cout<< "base "<<i->op(0)<<" degree "<<i->op(1)<<endl;
  }




 cout<<setw(48)<<"Overall dirivativ "<<latex<<(divpoly(dirac_trace(ine),dp,to_dif).simplify_indexed().subs(dp ==0)*indexed(dp,to_dif)).simplify_indexed(global_sp())<<endl;
  /////////////////////////////////////
  /// Convert to Rational expression///
  //ine = factor(ine).normal();//////////
  /////////////////////////////////////

  /*  cout<<"Rational function "<<ine<<endl;
  ex dn = ine.numer_denom();
  cout<<"Numer : "<< dirac_trace(dn.op(0)).expand()<<" \nand denom" <<dn.op(1)<<endl;
  /// diferentiation of Numerator
  cout<<setw(48)<<"Overall dirivativ "<<(divpoly(dirac_trace(dn.op(0)),dp,to_dif).simplify_indexed().subs(dp ==0)*indexed(dp,to_dif)).simplify_indexed()<<endl;
  /// Numereator wo matrices prepared for div
  ex womn = dirac_trace(dn.op(0));
  cout<<"\nRational div:"<< (((factor(dn.op(1).subs(dp==0))*factor(divpoly(womn,dp,to_dif).subs(dp==0)) - factor(womn.subs(dp==0))*factor(divpoly(dn.op(1),dp,to_dif).subs(dp==0))))*indexed(dp,to_dif)).simplify_indexed()<<endl;
  /// diferentiation of Denominator
//  divpoly(dn.op(1),dp);
*/

  //div(num_den[0]);
  //div(num_den[1]);

  cout << std::setw (10);
  cout << 77 << endl;
}

void afdcpp::mark_ext_vert(Graph& Gtmp )
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  Graph::vertex_iterator viu,viu_end;
  Graph::out_edge_iterator oei,oei_end;
  Graph G = Gtmp.parent();
  int n_ext = 0;
    cout<<"vertieces in parent :" <<num_vertices(Gtmp.root())<<"And in child:"<<num_vertices(Gtmp)<<endl;
    ///Set all vertices in internal state
  for (tie(viu, viu_end) = vertices(Gtmp);viu != viu_end; ++viu) 
          Gtmp[*viu].external = false;
  /// Marking vertieces as external for subgraph
  for (tie(viu, viu_end) = vertices(Gtmp);viu != viu_end; ++viu) 
    {
      if((out_degree(*viu,Gtmp) < out_degree( Gtmp.local_to_global(*viu), G))|| G[Gtmp.local_to_global(*viu)].external)
        {
          Gtmp[*viu].external = true;
          n_ext++;
          cout<<*viu<<" ####external found!!!!"<<endl;
          /// Vertex Momentum Flow Calculation
          /*         ex fp;
          for(tie(oei,oei_end) = out_edges(Gtmp.local_to_global(*viu),G);oei != oei_end;++oei)
            fp+=G[*oei].gmomentum;
          for(tie(oei,oei_end) = out_edges(vertex(*viu,Gtmp),Gtmp);oei!=oei_end;++oei)
            fp-=Gtmp[*oei].gmomentum;
            cout<<"  MOMENTUM: "<<fp<<endl;*/
          
        }
    }
  cout<<setw(20)<<std::left<<n_ext<<setw(30)<<"ext Vertices in graph"<<endl;
}
void afdcpp::M(Graph& Gtmp)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |                                     | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
/*
    Graph::vertex_iterator viu,viu_end;
  Graph::out_edge_iterator oei,oei_end;
  Graph G = Gtmp.parent();
  //  cout<<"vertieces in parent :" <<num_vertices(g.parent())<<"And in child:"<<num_vertices(g)<<endl;
  /// Marking vertieces as external for subgraph
      for (tie(viu, viu_end) = vertices(Gtmp);viu != viu_end; ++viu) 
        {
          if(out_degree(*viu,Gtmp) < out_degree( Gtmp.local_to_global(*viu), G) )
            {
              Gtmp[vertex(*viu,Gtmp)].external = true;
               
              cout<<"################external found!!!!"<<endl;
              /// Vertex Momentum Flow Calculation
              ex fp;
              for(tie(oei,oei_end) = out_edges(Gtmp.local_to_global(*viu),G);oei != oei_end;++oei)
                fp+=G[*oei].gmomentum;
              for(tie(oei,oei_end) = out_edges(vertex(*viu,Gtmp),Gtmp);oei!=oei_end;++oei)
                fp-=Gtmp[*oei].gmomentum;
              cout<<"  MOMENTUM: "<<fp<<endl;
              */
                ex e37 = ex_from_graph(Gtmp);
  symbol n("n");
/// Working example
 cout<< "\nExternal changed 1 " << e37.subs(k1-p1+wild() == n + wild())<<endl;
 
 
  cout<< "\nExternal changed 2 " << e37.subs(indexed(k1-p1 + wild(0),varidx(wild(1),4)) == indexed(n + wild(0),varidx(wild(1),4)))<<endl;
   cout<< "\nExternal changed 3 " << e37.subs(indexed(k1-p1 + wild(0),varidx(wild(1),4)) == indexed(n + wild(0),varidx(wild(1),4)))<<endl;
  

/*
              Gtmp[*viu].gmomentum = fp;
            }
        }
  */

  //ex mom = k1+p1;
  symbol k;
  ex mom = k;
  symbol sim = m0;

}
///  Get moments and fields of really external legs
///  vertex descriptor is global
ext_vert_mom_field_tuple_lst afdcpp::get_ext_legs(Graph& G)
{
  ext_vert_mom_field_tuple_lst l1;
  Graph::vertex_iterator viu,viu_end;
  Graph::vertex_descriptor vdtmp;
  for (tie(viu, viu_end) = vertices(G); viu != viu_end; ++viu)
    if(G[*viu].external == true)
      {

        std::list<pair<ex,field> >::iterator iter = ext_legs[make_pair(get_property(G.root(),graph_num),G.local_to_global(*viu))].begin();
        cout<<"Works get_ext "<<get_property(G.root(),graph_num)<<endl;
        std::list<pair<ex,field> >::iterator iter_end = ext_legs[make_pair(get_property(G.root(),graph_num),G.local_to_global(*viu))].end();
        for(iter;iter != iter_end;++iter)
          l1.push_back(make_tuple(G.local_to_global(*viu),iter->second,iter->first));
      }
  return l1;
}

// find momentums external for the diagram
//
// *** may be IN and OUT flags for momentums???
//
exvector afdcpp::get_ext_momentum(Graph& Gtmp)
{



  /// Vertex Momentum Flow Calculation
  Graph::vertex_iterator viu,viu_end,vi,vi_end;
  Graph::out_edge_iterator oei,oei_end;
  exvector e1;
  for(tie(vi,vi_end)=vertices(Gtmp);vi!=vi_end;++vi)
    if(is_external_vert(*vi,Gtmp))
      {
        e1.push_back(get_ext_mom_of_vertex(*vi,Gtmp));
        cout<<"Mom of vertex "<<Gtmp.local_to_global(*vi)<<" P= "<<e1.back()<<endl;
      }
  return e1;
  /* if(!Gtmp.is_root())
    {
  Graph G = Gtmp.root();
  ex fp;
  for(tie(oei,oei_end) = out_edges(Gtmp.local_to_global(*viu),G);oei != oei_end;++oei)
    fp+=G[*oei].gmomentum;
  cout<<"glob mom scaned"<<endl;
  for(tie(oei,oei_end) = out_edges(vertex(*viu,Gtmp),Gtmp);oei!=oei_end;++oei)
    fp-=Gtmp[*oei].gmomentum;
  cout<<"  MOMENTUM: "<<fp<<endl;
  
  ex e37 = ex_from_graph(Gtmp);
  symbol n("n");
  /// Working example
  cout<< "\nExternal changed 1 " << e37.subs(k1-p1+wild() == n + wild())<<endl;
  
  
  cout<< "\nExternal changed 2 " << e37.subs(indexed(k1-p1 + wild(0),varidx(wild(1),4)) == indexed(n + wild(0),varidx(wild(1),4)))<<endl;
  cout<< "\nExternal changed 3 " << e37.subs(indexed(k1-p1 + wild(0),varidx(wild(1),4)) == indexed(n + wild(0),varidx(wild(1),4)))<<endl;
  }*/
}

/// Get information about momentum of GLOBAL external vertex
std::pair<ex,bool> afdcpp::get_global_ext_mom_info(Graph::vertex_descriptor vdis,Graph& G)
{
  ex sum_mom = 0;
  if(G[vdis].external)
    {
       std::list<pair<ex,field> >::iterator iter = ext_legs[make_pair(get_property(G.root(),graph_num),G.local_to_global(vdis))].begin();
        cout<<"Works get_ext "<<get_property(G.root(),graph_num)<<endl;
        std::list<pair<ex,field> >::iterator iter_end = ext_legs[make_pair(get_property(G.root(),graph_num),G.local_to_global(vdis))].end();
        for(iter;iter != iter_end;++iter)
          sum_mom += (iter->first);
        return make_pair(sum_mom,true);
    }
  else
    {
      return make_pair(0,false);
    }
}

/// get summ of momentums external localy ol globaly to VDIS
ex afdcpp::get_ext_mom_of_vertex(Graph::vertex_descriptor vdis,Graph& Gtmp)
{
  Graph::vertex_iterator viu,viu_end,vi,vi_end;
  Graph::out_edge_iterator oei,oei_end;
  if(is_global_ext_vert(vdis,Gtmp))
    return get_global_ext_mom_info(vdis,Gtmp).first;
  if(is_local_ext_vert(vdis,Gtmp))
    {
      Graph G = Gtmp.root();
      ex fp;
      for(tie(oei,oei_end) = out_edges(Gtmp.local_to_global(vdis),G);oei != oei_end;++oei)
        fp+=G[*oei].gmomentum;
      cout<<"glob mom scaned"<<endl;
      for(tie(oei,oei_end) = out_edges(vertex(vdis,Gtmp),Gtmp);oei!=oei_end;++oei)
        fp-=Gtmp[*oei].gmomentum;
      cout<<"  MOMENTUM: "<<fp<<endl;
      return fp;
    }
}


/// SORTED list of fields on external lines, graph or SUBGRAPH
std::vector<field> afdcpp::get_ext_fields(Graph& G)
{
  if(G.is_root())
    {
      std::vector<field> l1;
      Graph::vertex_iterator viu,viu_end;
      Graph::vertex_descriptor vdtmp;
      for (tie(viu, viu_end) = vertices(G); viu != viu_end; ++viu)
        if(G[*viu].external == true)
          {
            std::list<pair<ex,field> >::iterator iter = ext_legs[make_pair(get_property(G,graph_num),*viu)].begin();
            std::list<pair<ex,field> >::iterator iter_end = ext_legs[make_pair(get_property(G,graph_num),*viu)].end();
            for(iter;iter != iter_end;++iter)
              l1.push_back(iter->second);
          }
      //cout<<"List of ext FIELDS constructed "<<l1<<endl;
      return l1;
    }
  else
    {
      std::vector<field> l1;
      Graph::vertex_iterator viu,viu_end;
      Graph::out_edge_iterator oei,oei_end;
      for (tie(viu, viu_end) = vertices(G); viu != viu_end; ++viu)
        if(is_external_vert(*viu,G))
          {
            // ext line lying in graph
            if(out_degree(*viu,G) < out_degree(G.local_to_global(*viu),G.root()))
              {
                //loop over rootgraph edges
                for(tie(oei,oei_end)=out_edges(G.local_to_global(*viu),G.root());oei!=oei_end;++oei)
                  //if(!find_edge(*oei,G))
                    {
                      if(source(*oei,G.root())==G.local_to_global(*viu))l1.push_back(G.root()[*oei].from_field);
                      if(target(*oei,G.root())==G.local_to_global(*viu))l1.push_back(G.root()[*oei].to_field);
                      }
              }
            // ext line of graph
            if(G.root()[G.local_to_global(*viu)].external)
              {
                std::list<pair<ex,field> >::iterator iter = ext_legs[make_pair(get_property(G.root(),graph_num),G.local_to_global(*viu))].begin();
                std::list<pair<ex,field> >::iterator iter_end = ext_legs[make_pair(get_property(G.root(),graph_num),G.local_to_global(*viu))].end();
                for(iter;iter != iter_end;++iter)
                  l1.push_back(iter->second);
              }
          }
      //cout<<"List of ext FIELDS constructed "<<l1<<endl;
      return l1;
    }
}
/// type,momentum,expr
tuple<int,ex,ex> afdcpp::T_expr_prop(field f1,field f2,int from,int to,ex mom)
{
  if(f1 == phi && f2 == phi)
    {
      prop<phi,phi> pr(from,to,mom);
      return make_tuple(pr.t,pr.momentum,pr.expr);
    }
}

//<momentum,expr>
std::pair<ex,ex> afdcpp::T_expr3(Graph::vertex_descriptor& vd,tuple_lst& tl,ext_vert_mom_field_tuple_lst& el)
{
  std::vector<field> f_lst; /// list of particles on vertex legs
  std::vector<ex> m_lst;    /// list of moments on vertex legs
  for(tuple_lst::iterator titer = tl.begin();titer != tl.end();++titer)
    {
      f_lst.push_back(get<2>(*titer));
      m_lst.push_back(get<3>(*titer));
    }
  for(ext_vert_mom_field_tuple_lst::iterator eiter = el.begin();eiter != el.end();++eiter)
    {
      f_lst.push_back(get<1>(*eiter));
      m_lst.push_back(get<2>(*eiter));
    }
  cout<<"MK: "<<f_lst.size()<<" "<<m_lst.size()<<endl;
  ///(phi,phi,phi,phi) vertex
  if(f_lst[0] == phi,f_lst[1] == phi,f_lst[2] == phi,f_lst[3] == phi)
    {
      vert4<phi,phi,phi,phi> v4(vd,m_lst[0],m_lst[1],m_lst[2],m_lst[3]);
      return make_pair(0,v4.expr);
    }
  else
    {
      cout<<"Unknown vertex=)"<<endl;
      symbol jkl;
      return make_pair(0,jkl);
    }
}

//<momentum,expr>
std::pair<ex,ex> afdcpp::T_expr4(Graph::vertex_descriptor& vd,tuple_lst& tl,ext_vert_mom_field_tuple_lst& el)
{
  std::vector<field> f_lst; /// list of particles on vertex legs
  std::vector<ex> m_lst;    /// list of moments on vertex legs
  for(tuple_lst::iterator titer = tl.begin();titer != tl.end();++titer)
    {
      f_lst.push_back(get<2>(*titer));
      m_lst.push_back(get<3>(*titer));
    }
  for(ext_vert_mom_field_tuple_lst::iterator eiter = el.begin();eiter != el.end();++eiter)
    {
      f_lst.push_back(get<1>(*eiter));
      m_lst.push_back(get<2>(*eiter));
    }
  cout<<"MK: "<<f_lst.size()<<" "<<m_lst.size()<<endl;
  ///(phi,phi,phi,phi) vertex
  if(f_lst[0] == phi,f_lst[1] == phi,f_lst[2] == phi,f_lst[3] == phi)
    {
      vert4<phi,phi,phi,phi> v4(vd,m_lst[0],m_lst[1],m_lst[2],m_lst[3]);
      return make_pair(0,v4.expr);
    }
  else
    {
      cout<<"Unknown vertex=)"<<endl;
      symbol jkl;
      return make_pair(0,jkl);
    }
}
/// pair<div_subraph,T(G)>
Graph afdcpp::K(Graph& Gdiv,Graph& GT)
{
  cout<<colors::blue<<"test in K"<<endl;
  if(!GT.is_root())cout<<colors::yelow<<"FAIL"<<endl;
  print(GT);
  cout<<colors::def_col;
  //Graph_1PI_lst l_div = make_subgraphs(GT);
  std::list<Graph> l_div;
  /*
  BOOST_FOREACH(Graph g,l_div)
    {
      if(isSubgraph(Gdiv,g))
        {
          cout<<"45 found div sub"<<endl;
          return g;
        }
    }
*/
  if(l_div.size())
    {
      for(std::list<Graph>::iterator it = l_div.begin();it!=l_div.end();++it)
        if(isSubgraph(Gdiv,*it))
          {
            cout<<"K found div sub"<<endl;
            //Graph g76 = *it;
            return *it;
          }
    }
  else
    {
      cout<<" MEGAERROR=)"<<endl;
      return 0;
    }
}


Graph afdcpp::T(Graph& Gtmp)
///////////////////////////////////////////////
///  _____________________________________  ///
/// |                                     | ///
/// |      Constructing T(G) terms        | ///
/// |_____________________________________| ///
///                                         ///
///////////////////////////////////////////////
{
  if(Gtmp.is_root())
    cout<<"\n ERROR: Can't apply T(G) on G, where G is a root graph object "<<endl;
  else
    {
      cout<<"and Gtmp"<<endl;
      print(Gtmp);
      Graph::vertex_descriptor gvdis,vdis;
      bool is_in_graph;
      Graph::vertex_iterator viu,viu_end;
      Graph::edge_iterator ei,ei_end;
      Graph::out_edge_iterator oei,oei_end;
      Graph G = Gtmp.root(); /// Main diagram graph

      cout<<"print Gtmp.root"<<endl;
      print(G);

      Graph GT;              /// Graph after T-operation
      get_property(GT,graph_num) = get_property(Gtmp,graph_num);
      bool from_is_in_div_graph,to_is_in_div_graph;
      bool inserted;
      Graph::edge_descriptor edis;
      // tuple<vdis,edis,field,momentum>
      tuple_lst tleg_lst;
      for (tie(ei, ei_end) = edges(G);ei != ei_end; ++ei)
        {
          //          tie(vdis,from_is_in_div_graph) = Gtmp.find_vertex(G[*ei].from_vertex);
          //          tie(vdis,to_is_in_div_graph) = Gtmp.find_vertex(G[*ei].to_vertex);

          from_is_in_div_graph =find_vert(G[*ei].from_vertex,Gtmp);
          to_is_in_div_graph = find_vert(G[*ei].to_vertex,Gtmp);
          
          if(from_is_in_div_graph && to_is_in_div_graph)                /// edge is in div subgraph
            {
                            cout<<"full edge in subgraph"<<endl;
            }
          if((!from_is_in_div_graph) && (!to_is_in_div_graph))            /// adding edge to T(G)
            {
              tie(edis, inserted) = add_edge(G[*ei].from_vertex,G[*ei].to_vertex,GT);
              GT[edis].type = G[*ei].type;
              GT[edis].gmomentum = G[*ei].gmomentum;
              GT[edis].from_vertex = G[*ei].from_vertex;
              GT[edis].to_vertex = G[*ei].to_vertex;
              GT[edis].gexpr = G[*ei].gexpr;
              cout<<"Inserted full edge"<<endl;
              
            }
          if((!from_is_in_div_graph) && to_is_in_div_graph)           /// from field     
            {
              tleg_lst.push_back(make_tuple(G[*ei].from_vertex,*ei,G[*ei].from_field,G[*ei].gmomentum));
              cout<<"To "<<G[*ei].to_vertex<<"is in div graph, and adding"<<G[*ei].from_vertex<<endl;
            }
          if(from_is_in_div_graph && (!to_is_in_div_graph))         /// to field     
            {
              tleg_lst.push_back(make_tuple(G[*ei].to_vertex,*ei,G[*ei].to_field,G[*ei].gmomentum));
              cout<<"From "<<G[*ei].from_vertex<<"is in div graph, and adding"<<G[*ei].to_vertex<<endl;
            }
        }
      /// Adding vertex information
      /// Check if no vertices in GT
      ///
      /// !!!!!!!!!!SINGLE VERTICES ARE NOT ADDED!!!!!!!!!!!!
      ///
      /// Adding info in all GT vertices
      for (tie(viu, viu_end) = vertices(GT);viu != viu_end; ++viu)
        {
          gvdis = vertex(*viu,G); // vertex discriptor of root graph vert
          GT[gvdis].gexpr = G[gvdis].gexpr ;
          GT[gvdis].gmomentum = G[gvdis].gmomentum ;
          GT[gvdis].external = G[gvdis].external;
          cout<<"GT has "<<num_vertices(GT)<<" vertices\n";
          
        }
      ext_vert_mom_field_tuple_lst l3 = get_ext_legs(Gtmp);
      cout<<"FNumber of added verts: "<<tleg_lst.size()<<" and ext in GT "<<num_ext_vert(GT)<<" and real ext in M = "<<l3.size()<<" verts to shrink: "<< l3.size()+tleg_lst.size()<<endl;
      /// number of external legs
      int n_ext_l = l3.size();
      /// number of internal legs
      int n_int_l = tleg_lst.size(); 
      /// number of legs in subgraph to be shrinked
      int n_legs = n_ext_l+n_int_l;
      switch(n_legs)
        {
        case 2: // ins propagator
          {
            cout<<"\n\n  ext legs: "<<n_ext_l<<"       I want propagator!!!\n\n"<<endl;
            
          }
        case 3: // ins 3 part vertex
          {
             cout<<"\n\n  ext legs: "<<n_ext_l<<"       I want 3 part vert!!!\n\n"<<endl;
          }
        case 4: // ins 4 part vertex
          {
            vdis = Gtmp.local_to_global(vertex(0,Gtmp));
            cout<<" VDIS:"<<vdis<<endl;
            if(n_int_l > 0)   /// if edges between new vertex and old exists
              for(tuple_lst::iterator titer = tleg_lst.begin();
                  titer != tleg_lst.end();++titer)
                {
                  cout<<"Before ins"<<endl;
                  print_edges2(GT, get(vertex_index, GT), get(edge_index, GT));
                  print_vertices(GT,get(vertex_index, GT));
                  //print_edges(GT,get(edge_index, GT));
                  print_graph(GT);
                  /// add edge from new vertex (vdis) to vertex <0> in half of edge
                  tie(edis, inserted) = add_edge(vdis,get<0>(*titer),GT);
                  if(get<0>(*titer) == G[get<1>(*titer)].from_vertex)//edge from M(G)
                    {
                      tie(GT[edis].type,GT[edis].gmomentum,GT[edis].gexpr)=
                        T_expr_prop(G[get<1>(*titer)].from_field,
                                    G[get<1>(*titer)].to_field,
                                    G[get<1>(*titer)].to_vertex,vdis,
                                    G[get<1>(*titer)].gmomentum);
                      GT[edis].from_vertex = get<0>(*titer);
                      GT[edis].to_vertex = vdis;
                      GT[edis].from_field = G[get<1>(*titer)].from_field;
                      GT[edis].to_field = G[get<1>(*titer)].to_field;
                      /// adding expression to vertex from root graph
                      GT[get<0>(*titer)].gexpr = G[get<0>(*titer)].gexpr;
                      GT[get<0>(*titer)].external = G[get<0>(*titer)].external;
                      GT[get<0>(*titer)].gmomentum = G[get<0>(*titer)].gmomentum;
                    }
                  if(get<0>(*titer) == G[get<1>(*titer)].to_vertex) // edge to M(G)
                    {
                      tie(GT[edis].type,GT[edis].gmomentum,GT[edis].gexpr)=
                        T_expr_prop(G[get<1>(*titer)].from_field,
                                    G[get<1>(*titer)].to_field, vdis,
                                    G[get<1>(*titer)].to_vertex,
                                    G[get<1>(*titer)].gmomentum);
                      GT[edis].from_vertex = vdis;
                      GT[edis].to_vertex = get<0>(*titer);
                      GT[edis].from_field = G[get<1>(*titer)].from_field;
                      GT[edis].to_field = G[get<1>(*titer)].to_field;
                      /// adding expression to vertex from root graph
                      GT[get<0>(*titer)].gexpr = G[get<0>(*titer)].gexpr;
                      GT[get<0>(*titer)].external = G[get<0>(*titer)].external;
                      GT[get<0>(*titer)].gmomentum = G[get<0>(*titer)].gmomentum;
                    }
                  
                  cout<<"After ins"<<endl;
                  print_vertices(GT,get(vertex_index, GT));
                  //print_edges(GT,get(edge_index, GT));
                  print_graph(GT);
                  cout<<"-----"<<endl;
                }
            else        /// new edges are no added diagram==single vertex
              {
               vdis =  add_vertex(GT);
               cout<<"Single vert "<<endl;
                  print_graph(GT);
                  cout<<"-----"<<endl;
              }
            
            // if external legs exists in, than it is external vertex
            GT[vdis].external = (n_ext_l > 0);
            tie(GT[vdis].gmomentum,GT[vdis].gexpr)=T_expr4(vdis,tleg_lst,l3);
            //GT[get<0>(*titer)].gmomentum = G[get<0>(*titer)].gmomentum;
            cout<<ex_from_graph(GT)<<endl;

            /*tie(edis, inserted) = add_edge(G[*ei].from_vertex,G[*ei].to_vertex,GT);
              GT[edis].type = G[*ei].type;
              GT[edis].gmomentum = G[*ei].gmomentum;
              GT[edis].from_vertex = G[*ei].from_vertex;
              GT[edis].to_vertex = G[*ei].to_vertex;
              GT[edis].gexpr = G[*ei].gexpr;*/
          }
        }
      cout<<"print T(G)"<<endl;
      print(GT);
        
      cout<<"Creates T(G)"<<endl;
      return GT;
    }
      
            

/*
  //  cout<<"vertieces in parent :" <<num_vertices(g.parent())<<"And in child:"<<num_vertices(g)<<endl;
  /// Marking vertieces as external for subgraph
      for (tie(viu, viu_end) = vertices(Gtmp);viu != viu_end; ++viu) 
        {
          if(out_degree(*viu,Gtmp) < out_degree( Gtmp.local_to_global(*viu), G) )
            {
              Gtmp[vertex(*viu,Gtmp)].external = true;
               
              cout<<"################external found!!!!"<<endl;
              /// Vertex Momentum Flow Calculation
              ex fp;
              for(tie(oei,oei_end) = out_edges(Gtmp.local_to_global(*viu),G);oei != oei_end;++oei)
                fp+=G[*oei].gmomentum;
              for(tie(oei,oei_end) = out_edges(vertex(*viu,Gtmp),Gtmp);oei!=oei_end;++oei)
                fp-=Gtmp[*oei].gmomentum;
              cout<<"  MOMENTUM: "<<fp<<endl;
              */
                ex e37 = ex_from_graph(Gtmp);
  symbol n("n");
/// Working example
 cout<< "\nExternal changed 1 " << e37.subs(k1-p1+wild() == n + wild())<<endl;
 
 
  cout<< "\nExternal changed 2 " << e37.subs(indexed(k1-p1 + wild(0),varidx(wild(1),4)) == indexed(n + wild(0),varidx(wild(1),4)))<<endl;
   cout<< "\nExternal changed 3 " << e37.subs(indexed(k1-p1 + wild(0),varidx(wild(1),4)) == indexed(n + wild(0),varidx(wild(1),4)))<<endl;
  

/*
              Gtmp[*viu].gmomentum = fp;
            }
        }
  */

  //ex mom = k1+p1;
  symbol k;
  ex mom = k;
  symbol sim = m0;

}

/// Generation of all 1PI subgraphs of graph






